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Abstract 



In this paper we develop a family of three 8-step methods, optimized for 
the numerical integration of oscillatory ordinary differential equations. We 
have nullified the phase-lag of the methods and the first r derivatives, where 
r = {1,2,3}. We show that with this new technique, the method gains 
K^ \ efficiency with each derivative of the phase-lag nullified. This is the case for 

^ ■ the integration of both the Schrodinger equation and the N-body problem. 

Q\ • A local truncation error analysis is performed, which, for the case of the 

^ . Schrodinger equation, also shows the connection of the error and the energy, 

revealing the importance of the zero phase-lag derivatives. Also the stability 

00 \ analysis shows that the methods with more derivatives vanished, have a 

^ ' bigger interval of periodicity. 
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1. Introduction 



The numerical integration of systems of ordinary differential equations 
with oscillatory solutions has been the subject of research during the past 
decades. This type of ODEs is often met in real problems, like the N-body 
problem and the Schrodinger equation. 

There are some special techniques for optimizing numerical methods. 
Trigonometrical fitting and phase-fitting are some of them, producing meth- 
ods with variable coefficients, which depend on f = uh, where u is the 
dominant frequency of the problem and h is the step length of integration. 

For example Raptis and Allison have developed a two-step exponentially- 
fitted method of order four in 19| and Kalogiratou and Simos have con- 



structed a two-step P-stable exponentially-fitted method of order four in 



13l |. Also Panopoulos, Anastassi and Simos have constructed two optimized 



16 



eight-step methods with high or infinite order of phase-lag in 

Some other notable multistep methods for the numerical solution of oscil- 
lating IVPs have been developed by Chawla and Rao in [G], who produced a 
three-stage, two-Step P-stable method with minimal phase-lag and order six 
and by Henrici in |9| , who produced a four-step symmetric method of order 
six. Also some recent research work in numerical methods can be found in 

Q, B, i, 0, H, 0, 0> B, B, 0, H, and [l7 . 

Trigonometrically fitted methods of high trigonometric order are well 
known for their high efficiency in the integration of the Schrodinger equation, 
especially when using a high value of energy. However higher trigonometric 
order is not rendering them more efficient for all types of oscillatory prob- 
lems. On the other hand, phase-lag does not give us the opportunity to 
provide such methods, that for example perform well when integrating the 
Schrodinger equation for high values of energy. 

In this paper we present a methodology for optimizing numerical meth- 
ods, through the use of phase-lag and its derivatives with respect to v. 
More specifically, given a classical (i.e. with constant coefficients) numer- 
ical method, we can provide a family of optimized methods, each of which 
has zero {PL} or zero {PL and PL'} or zero {PL, PL' and PL"} etc. 

With this new technique we provide methods that perform well during 
the integration of the Schrodinger equation for high values of energy, but also 
that perform well on other real problems with oscillatory solution, like the 
N-body problem. 



2. Phase-lag and stability analysis of symmetric multistep methods 

For the numerical solution of the initial value problem 

y" = f{x,y) (1) 

multistep methods of the form 

•m m 

^ aiVn+i = h^^ bif{Xn+i, Vn+i) (2) 

j=0 i=0 

with m steps can be used over the equally spaced intervals {xj}™ q G [a, h] 
and h = |xj+i — Xi\, i = 0(l)m — 1. 

If the method is symmetric then Oj = a^-i and 6j = bm-i, i = 0(1) [yj. 
Method ([2]) is associated with the operator 

m m 

L{x) = 2_. ciiu{x + ih) — h^ 2_. biu"{x + ih) (3) 

i=0 1=0 

where u E C"^. 

Definition 1. The multistep method ([3]) is called algebraic of order p if 
the associated linear operator L vanishes for any linear combination of the 
linearly independent functions 1, x, x^, . . . , x^^^. 

When a symmetric 2A;-step method, that is for i = —k{l)k, is applied to 
the scalar test equation 

y" = -u^'y (4) 

a difference equation of the form 

Ak{v)yn+k + ••• + Ai{v)yn+i + Ao{v)yn 

+Ai{v)yn-i + ... + Akiv)yn-k = (5) 

is obtained, where v = uh, h is the step length and Ao{v), Ai{v), . . ., Ak{v) 
are polynomials of v. 

The characteristic equation associated with ([5]) is 

Ak{v)s'' + ... + Ai{v)s + Ao{v) + Ai{v)s-^ + ... + Ak{v)s^'' = (6) 



Theorem 1. 20| The symmetric 2k-step m,ethod with characteristic equa- 



tion given by (jS]) has phase-lag order q and phase-lag constant c given by 



,+2,r,( 1+^) = 2Afc(^) <^o<kv) + ... + 2Aj{v) cos(jt;) + ... + Ao{v) 

^ ' 2k^Ak{v) + ... + 2fA,{v) + ... + 2A,{v) ^' 

The formula proposed from the above theorem gives us a direct method 
to calculate the phase-lag of any symmetric 2k- step method. 

The characteristic equation has m characteristic roots Aj, i = 0(l)m — 1. 

Definition 2. [14!] If the characteristic roots satisfy the conditions |Aj| ^ 
1,2 = 0{l)m—l for all s = 6h, then we say that the method is unconditionally 
stable. 

Definition 3. [14] If the characteristic roots satisfy the conditions Ai = 
e^<^W, As = e-^-^^") |Ai| < 1, i = 3(l)m - 1 for all s < sq, where s = Oh 
and (f){s) is a real function of s, then we say that the method has interval of 
periodicity (0, s^). 



Definition 4. [1^ Method ([2]) is called P-stable if its interval of periodicity 
is (0, oo). 

3. Construction of the new optimized multistep methods 



We consider the multistep symmetric method of Quinlan-Tremaine 18 
with eight steps and eighth algebraic order: 



(8) 



2/4 = -y-A - ctsivs + y-s) - ^2(2/2 + y~2) - ai{yi + y~i) 

+ h' (63(/3 + /-3) + &2(/2 + /-2) + 6l(/l + /-l) + 60/0) 

where 

0-3 = —2, ^2 = 2, ai = —1, 

17671 , 23622 , 61449 , 50516 ,^. 

Oq = , Oo = , 61 = , On = , (9) 

^ 12096' 12096' 12096' 12096 

yi = y{x + ih) and fi = f{x + ih, y{x + ih)) 

We also consider the optimized method, that is based on the above one, 
with zero phase-lag constructed by Panopoulos, Anastassi and Simos in [16| . 
The coefficients are given below: 



on. 601 ^ .,109 , ,,, 101 

bo = -20b3 + —-, b2 = -6h + --, 6i = 15 63--r 
24 16 6 

63 = , where 

9Q D 

C = -192 (cos (v))'^ + 192 (cos (t;))^ + (96 - 327t;2) (cos {v)f 
+ (-120 + 404 1;2) cos (v) - 137 1;^ + 24 
D = v^ {cos (v) - if , 

where v = uh and the a, coefficients remain the same. The Taylor series 
expansions of the coefficients are: 



12629 45767 o 164627 . 520367 ^ 



3024 36288 2395008 792529920 
76873 „ 9190171 ^ 6662921 ,^ 



4483454976 160059342643200 1703031405723648 

_ 20483 45767 2 164627 4 520367 g 

^ ~ 4032 ~ 48384 ^ ^ 3193344 ^ ~ 1056706560 ^ 

76873 g 9190171 ^q 6662921 ^^ 

^5977939968^ ~ 213412456857600 ^ ~ 2270708540964864 ^ 



3937 45767 ^ 164627 , 520367 
^2 ^ 



62 = + v^ - V* + —— v^ 



2016 120960 7983360 2641766400 
76873 g 9190171 ^q 6662921 ^^ 



14944849920 ^ 533531142144000 ^ 5676771352412160 ^ 



17671 45767 ^ 164627 . 520367 
^3 



63 = i::::^, — -■ ^2^ ^4 ^e 

12096 725760 47900160 15850598400 

76873 g 9190171 ^q 6662921 ^^ 

^89669099520^ ~ 3201186852864000^ ~ 34060628114472960^ 

We want to produce three new methods that, apart from zero phase-lag, 
will also have zero r derivatives of the phase-lag, where r = {1,2,3}. In 
particular the three new methods must satisfy these equations: 

• First method: {PL = 0, PL' = 0} 

• Second method: {PL = 0, PL' = 0, PL" = 0} 



• Third method: {PL = 0, PL' = 0, PL" = 0, PL'" = 0} 

Since we have four free coefficients 6j, i = {0, 1, 2, 3} (oj remain the same), 
the rest of the coefficients for each method will be determined by the algebraic 
conditions. 

3.1. First optimized method with zero PL and PL' 

The first method must satisfy the conditions {PL = 0,PL' = 0}, thus 
we need two coefficients to be determined by the maximum algebraic order. 

We use formula ([7j) to compute the phase-lag and then its first derivative 
in respect to v. 
where v = uj h, cu is the frequency and h is the step length used. 

PL = (96 (cos {v)f + (-96 + 48 v%3) (cos {v)f + 

(-48 + 24 v%2) (cos {v)f + (60 + (125 - 144 63 - 48 62) v^) cos (f ) 

-12 + (24 62 -95 + 96 63) v^) 
/(60 + 125t;2) 



PL' = -(- 4800 V (cos (v))* + (1152 631; - 9600 sin (v) v^ + 4800 v 
-4608 sin (v) ) (cos {v)f + I -3600 ( — + ^^ ) (-2 + ^^^^3) sin (v) 

+ (576 62 + 2400) v) (cos {v)f + ( -1200 ( — + vA {v%2 - 2) sin {v) 

-3456 17(63 + 1/3 62)) cos (1;) 

+3600 (g + v^) (-1 + (-^ + ^3 + ^ h) v^) sin (.) 

+2304 1; (-§ + 63 + 1/462) ) (12 + 25 v^)'^ 

The four equations to be solved are: 

95 125 

PL = 0, PL' = 0, 60 = -— + I663 + 662, 61 = — --963-462 

6 12 

and the coefficients are given below: 



OQ,num 


61 = 


48L> ' 


62 = 


"2,num 


&3 = 


"3,n«m 


12D ' 


24D ' 


481) 



&0 

where D = ((cos (v)) — 2 (cos (f )) + 2 cos (v) — l) f ^ and 

6 



(10) 



bo,num = —288 (cos (f )) V + 576 sin (v) (cos (f )) + 192 (cos (f )) v 
-192 sin (t;) (cos (v))^ + 190 (cos {v)Y v^ + 720 (cos (v))^ v 
-120 (cos (t;))^ V - 672 sin (w) (cos {v)f + 370 (cos {v)f v^ 
+168 sin (t;) (cos {v)f - 540 (cos {v)f v + 145 (cos {v)f v^ 
+168 cos (w) sin (v) — 70 cos (f ) f ^ — 72 cos (f ) f + 108 f 
-48 siniv) - 35 w^ 



bi,num = —768 (cos (f )) f + 1536 sin (v) (cos (f )) + 192 (cos (f )) v 

+500 (cos (t;))^ t;=' - 192 sin (v) (cos {v)f + 2400 (cos (w))^ v 

+ 1000 (cos (i;))^ t;3 - 2112 sin {v) (cos {v)f - 1980 (cos (t;))^ v 

+595 (cos (v))^ t;3 + 288 sin (v) (cos {v)f - 100 cos (v) t;^ 

+648 cos (t;) sin (v) - 192 cos (w) t; + 348 1; - 195 w^ - 168 sin (v) 



b2,num = —96 (cos (f )) f + 192 sin {v) (cos (f )) — 192 (cos (v)) v 
+ 192 sin (v) (cos (t;))^ + 624 (cos (v))^ v + 216 (cos {v)f v 
-480 sin (v) (cos (t;)f + 250 (cos {v)f v^ - 612 (cos {v)f v 
+215 (cos (t>)) t>^ — 72 sin (t>) (cos (f)) — 70 cos (f ) f ^ 
+216 cos (v) sin (v) — 24 cos (v)v + 8Av — 48 sin (v) — 35 f^ 



bs nura = -rz- 192 (cos {v)f V + 192 sin {v) (cos (t;))^ + 288 (cos {v)f v 

3 3 2 

+192 (cos (f)) t> — 192 sin (t>) (cos (f)) — 96 sin (t>) (cos (f )) 

-324 (cos {v)f V + 125 (cos {v)f v^ + 120 cos {v) sin (w) + 60 cos {v) v^ 

-24 sin(t;) + 36t;-65t;3 

The Taylor series expansions, used when f ^ 0, are given below: 



12629 45767 . 11483491 . 112258001 . 

fen = \ V^ V^ -\ V^ 

3024 18144 23950080 2615348736 
1703481341 „ 5614773343 ^ 10940565121 -,. 



784604620800 80029671321600 6307523724902400 

20483 45767 , 10476617 , 45578707 . 

4032 24192 31933440 1585059840 
1514526707 § 5016343559 ^q 19742264573 ^^ 

^1046139494400^ ~ 106706228428800 ^ ^ 17466988776652800^ 

3937 45767 o 1491199 . 321593093 . 

62 = \ v^ v^ H v^ 

2016 60480 15966720 43589145600 

189532561 o 460150601 ^ 28082396599 .^ 



523069747200 38109367296000 113535427048243200 



17671 45767 o 96865 , 21971953 

'J3 



12096 362880 19160064 261534873600 

82561 8 17608099 ^p _ 1184824691 ^2 _ 

+ 448345497600^ ~ 123122571264000 ^ ^ 75690284698828800 ^' 

^.^ Second optimized method with zero PL, PL' and PL" 

The second method must satisfy the conditions {PL = 0, PL' = 0, PL" = 
0}, thus we need one coefficient to be determined by the maximum algebraic 
order. 

We use formula ([7]) to compute the phase-lag and then its ffist and second 
derivative in respect to v. 



PL = (16 (cos (v))^ + (8 bav^ - 16) (cos {v)f + (4 biV^ - 8) (cos {v)y 
+ (10 + (-6 63 + 2 61) v^) cos (f ) - 2 + (-4 62 - 2 61 - 2 63 + 5) v^) 
/{l0 + {18b3 + 8b2 + 2bi)v^) 



^3 

V 



PL' = (-16 1; (9 63 + 4 62 + h) (cos {v)f + ((-160 + (-32 h^ - 128 62 
-288 63) v"^) sin iv) + 16 (f 63 + 4 62 + h) v) (cos (t;))^ + 

(-12 (5 + (9 63 + 4 62 + &i) t^') (-2 + M^) sin (^) 

+8 (9 63 + f &2 + &i) t^) (cos {v)f + (-4 (62t;' - 2) (5+ 
(9 63 + 4 62 + &i) t'^) sin iy) - 40 1; (3 63 + 62)) cos {y) - 

(5 + (-3 63 + fei) ^') (5 + (9 63 + 4 62 + h) v^) sin (t;) 

-8t;(|62-&3 + fei-f)) 

/(5 + (9 63 + 4 62 + &i)t^')' 



PL" = W9 ((-10368 (63 + I 62 + I 6i)^^ 

+3888 (i 61 - I + 63 + I 62) (63 + I &2 + I h) v^ - 3200 - 720 h. 

-320 62 - 80 61) (cos (^;))^ + (10368 (| + (63 + | &2 + | &i) t^') 

(fea + I ^2 + I &i) sin (t;) - 2916 (63 + | 62 + | 61)' 63^^^ 

+2592 (63 + I 62 + I 61) (62 + &3 + i &i) t^' + (-768 62' 

+ (2880 - 3936 63 - 384 61) 62 - 4968 63^ + (5580 - 984 61) 63 

+720 61 - 48 61^) t;2 + 1800 + 920 63 + 320 62 + 80 61) (cos {v)f 

+ (-9936 (&3 + i &i + i &2) (I + (&3 + I &2 + I &i) v"") V sin (t;) 

-648 62 (&3 + I &2 + I &i)' t;' + 9072 (§62 + 63 + 1 61) 

(63 + I 62 + I 61) ^^' + (-624 62' + (4280 - 2268 63 - 252 61) 62 

-1944 (63 + i6i) (63 + |6i-^))t^' + 2800 + 260 62 + 40 6i 

+360 63) (cos (t;))' + (-2592 (63 + i 61 + if 62) (§ 

+ (^3 + §62 + |6i) t;2) t;sin(t;) + 2187 (63 - ^ 61) (63 + | 62 + |6i)' 

-1863 (H 62 + i 61 + 63) (63 + I 62 + I 61) v^ + (480 62' 

+ (-2120 + 120 61 + 2520 63) 62 + 3240 63^ + (-4095 + 360 61) 63 

-555 61) t;2 - 1325 - 600 63 - 200 62) cos (t;) + 2160 (| + (63 + | 62 

+ |6i) t;2) t; (I62 + 63) sin(t;) + 324^2 {h + | 62 + |6i)'t;' 

-648 (-1 62 + 63 + I 61) (63 + I 62 + I 61) t^" + (144 62' 

+ (-520 + 22863 + 13261)62-216 (63 + I 61) (63-6i + ^))t;2 

-75 - 60 62 - 40 61 + 40 63) 

/(i+(^3 + |62 + |60^;^)^ 



\3 



,3 



The four equations to be solved are: 

PL = 0, PL' = 0, PL" = 0, 60 = 5 - 2 62 - 2 61 - 2 63 



and the coefficients are given below: 

7 "0,nMm , "l,rtM?Tt r ^2,num t ^3,num 

^°=^^' ^' = '^^' ^' = ^^' ^' = '^^' (11) 

where D = v"^ (sin (v)) (cos (v) — 1) and 

bo,num = -6 + 25 (cos {v)f v^ + 16 (cos {v)y v^ - 120 (cos (v))^ 

— 32 sin (t>) t> (cos (f)) — 96 sin (t>) t> (cos (f )) + 32 (cos (f )) v^ 

—36 cos (f ) f ^ + 15 cos (f ) f*^ + 20 (cos (f )) f "^ — 96 (cos (f )) 

+30 v^ (cos {v)f + 20 sin (v) t; - 12 1;^ + 10 (cos {v)f v^ 

+160 sin (v) V (cos (v)) + 140 sin (f ) v (cos (f )) 

— 60 sin (f) t> (cos (f)) — 134 sin (f) f (cos (f )) + 2 sin (f ) f cos (f ) 

+ 18 cos (v) + 30 (cos {v)f - 54 (cos (v))^ + 192 (cos {v)f 

+36 (cos (t;))^ + 24 (cos {v)f v^ - 64 (cos {v)f v^ + 88 (cos (v))^ v^ 

-68 (cos(t;)f t;2 + 20 (cos(t;))%2 



&i,n^.m = -18 - 192 (cos iv)Y + 120 (cos (t;))^ v'^ + 128 (cos (t;))^ t;^ 

-480 (cos {v)f - 320 sin (v) v (cos (t;))^ - 192 sin (v) v (cos {v)Y 

+64 (cos (t;))^ v^ - 104 cos (t;) v^ + 30 cos (t;) v^ + 15 1;^ 

+60 (cos (f)) f"^ — 192 (cos (w)) + 75 w"^ (cos (f)) + 64 sin (f ) t; 

—40 v^ + 496 sin (t>) f (cos (f )) + 680 sin {v) v (cos (w)) 

—320 sin (w) v (cos (f)) — 418 sin (f ) f (cos (v)) + 10 sin (f) f cos (f) 

+42 cos (v) + 162 (cos {v)f - 258 (cos {v)f + 528 (cos {v)f 

+408 (cos (t;))^ + 32 (cos {v)f v^ - 176 (cos {v)f v^ + 336 (cos (v))^ v^ 

-360 (cos {v)f v^ + 120 (cos (v))^ v^ 



b2,num = -6 - 96 (cos {v)Y + 15 (cos {v)f V* + 48 (cos {v)Yv^ 

—84 (cos (f )) — 128 sin (v) v (cos (f )) — 40 cos {v) v"^ + 15 cos {v) v^ 

+30 v^ (cos {v)f + 20 sin (t;) t; - 8 1;^ + 48 sin {v) v (cos {v)f 

+240 sin (f ) t> (cos (y)) — 48 sin (t>) f (cos (f )) 

—126 sin iy) v (cos (f )) — 6 sin (v) v cos (f ) + 18 cos (f ) + 42 (cos (f )) 

-114 (cos {v)Y + 48 (cos {v)f + 192 (cos {v)Y + 16 (cos (t;))^ v^ 

+ 128 (cos (t;)f t;2 - 136 (cos {v)Y v^ - 8 (cos {v)Y v^ 



10 



b3,num = 48 (cos {v)fv^ - 48 (cos {v)f + 48 (cos {v)f 

— 80 sin (t>) t> (cos (f)) — 48(cos(f)) f^ + 80 sin (t>) i; (cos (f)) 

-96 (cos (v))^ v^ + 72 (cos (v))^ + 96 (cos {v)f v^ - 78 (cos {v)f 

+104 sin {v) V (cos {v)) — 18 (cos (f )) — 102 sin (v) v (cos (f )) 

+48 (cos (f )) f ^ + 5 w^ (cos (f )) — 18 sin (f ) v cos (f ) — 48 cos (v) t>^ 

+30 cos (f ) + 10 cos (f ) f "^ — 6 + 16 sin {v)v + 5 v'^ 

The Taylor series expansions of the coefficients are given below: 



12629 45767 o 9837221 , 153204313 ^ 
" 3024 12096 7983360 653837184 

2356782689 „ 20347993339 ^ 8744186458121 -,. 



87178291200 9700566220800 77410518441984000 

20483 45767 ^ 2943449 , 107557349 ^ 

hi = V -\ V -u 

4032 16128 3548160 792529920 

5074066909 §_ 10190684747 ^q 5994017812967 ^^ _ 
^348713164800^ ~ 9484998082560 ^ ^ 103214024589312000^ 

3937 45767 o 8607 . 51408821 « 

H V ; V + 



2016 40320 39424 2724321600 

35318011 o 3348191339 ^ 56104711163 ,. 



34871316480 118562476032000 43667471941632000 

_ 17671 45767 2 22153 4 _ 41092123 g 

' ~ 12096 ~ 241920 ^ ^ 4561920 ^ ~ 130767436800 ^ 

7321421 „ 5642643317 ^ 210863655707 12 



348713164800 2134124568576000 681212562289459200 

3.3. Third optimized method with zero PL, PL' , PL" and PL'" 

All four free coefficients of the third method will be determined by con- 
ditions {PL = 0, PL' = 0, PL" = 0, PL'" = 0}. 

We use formula (jTj) to compute the phase-lag and then its ffist, second 
and third derivative in respect to v. 

PL = (16 (cos (v))^ + (-16 + 8 bsv"^) (cos {v)f + (4 63^^^ - 8) (cos {v)f 
+ (10 + (2 61 - 6 63) v^) cos (t;) - 2 + (-2 62 + &o) v^) / 
(10 + (18 63 + 8 62 + 2 61) t;2) 

11 



PL' = (-16 i; (9 63 + 462 + h) (cos (v))^ + ((-160 + (-128 62 
-32 61 - 288 63) v"^) sin (t;) + 16 (4 62 + 61 + f 63) v) (cos (t;))^ + 
(-12 (-2 + 63^') (5 + (9 63 + 4 62 + &i) t^') sin (t;) 
+8 (f 62 + 9 63 + &i) ^) (cos iv)f + (-4 (5 + (9 63 + 4 62 + &i) v^) 
(-2 + M') sin (v) - 120 1^ (63 + | &2)) cos (v) - (5 + (9 63 + 4 62 + &i) ^') 
(5 + ih - 363)^') sin(t;) + 5v{bo + f 63 - | 62 + f&i)) 
(5 + (9 63 + 4 62 + fei)t^')' 



PL" = ^ ((-2048 (|63 + &2 + l/46i)%^ + 768 (i 61 - f + 62 + | ^3) 

(f &3 + &2 + i &i) v^ - 320 62 - 80 61 - 3200 - 720 63) (cos (v))^ 
+ (2048 (I + (I 63 + 62 + i 61) v^) v{lh + h + l 61) sin (.;) 
-57663 (1^3 + b2 + lhyv^ + 1152 (1^3 + 62 + \h) (63 + 62 + i^i) v^ 
+ (-768 62^ + (2880 - 384 h - 3936 63) 62 - 4968 63^ + (5580 - 984 61) 63 
-48 61^ + 720 61) v^ + 1800 + 920 63 + 80 61 + 320 62) (cos {v)f 
+ (-1536 (f + (1 63 + &2 + i 61) ^') (i &i + f &3 + ^2) V sin (t;) 
-12862 (I63 + b2 + i6i)%^ + 1472 (§63 + X5^ + ^2) 
(i 63 + 62 + i 61) ^' + (-624 62' + (-252 61 + 4280 - 2268 63) 62 
-1944 (63 + I 61) (63 + 1 61 - ^)) t;2 + 2800 + 360 63 + 40 61 + 260 62) 
(cos {v)r + (-832 (I + (1 63 + 62 + 1 61) v^) (62 + g 63 + ^ 61) V sin (v) 
+432 (63 - ^61) (I63 + 62 + \h)\^ - 848 (If 63 + 62 + ^&i) 
(f 63 + 62 + i 61) v^ + (480 62' + (-2120 + 2520 63 + 120 61) 62 + 3240 63' 
+ (360 61 - 4095) 63 - 555 61) v^ - 1325 - 200 62 - 600 63) cos (v) 
+320 (62 + 3 63) (I + (I 63 + 62 + i 61) v^) V sin (v) + 64 62 
(1^3 + 62 + lhyv' + 32 (I63 + 62 + lh) (-963 + 62 - 61) v^ 
+ (24 62^ + (-220 - 18 61 - 162 63 - 60 60) 62 - 486 (63 + | 61) 
(&3 + ^ 60 + i + l^i)) v^ - 200 + 10 61 - IO62 + 9O63 + 25 60) 
/(|+(!63 + 62 + i6i)t;^)' 
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PL'" = (768 V ((9 63 + 4 ^2 + hf t;^ - i (9 63 + 4 62 + h) 
(61 - 40 + 9 63 + 4 62) ^' + 5 62 + f &i + f &3 + 25) 
(9 63 + 4 62 + h) (cos (t;))^ + (512 (5 + (9 63 + 4 62 + ^i) v^) 
((963 + 462 + Wfv^ - I (963 + 462 + 61) (61 - f + 963 + 462) t;2 



63 + f 62 + 25 + f 61) sin {v) - 432 ((9 63 + 4 62 + 61 



V 



-| (&i-f + 963 + 462)(963 + 462 + 6i)t;2 + f 62 + f 61 + 2O63 + 25) 
v (4 62 + 61 + f 63)) (cos {v)f + (108 (63 (9 63 + 4 62 + 61)' t;6 - 2 (4 63 
+4 62 + 61) (9 63 + 4 62 + 61) t;4 + (414 63' + (328 62 - 155 + 82 61) 
63 + 4 (61 + 4 62) (61 - 5 + 4 62)) t;2 - 50 - f 62 - f 61 - ^ 63) 
(5 + (9 63 + 4 62 + 61) v^) sin iv) - 672 v ((9 63 + 4 62 + 61)' 
(61 + fi 62 + 9 63) t;4 - i (9 63 + 4 62 + 61) (81 63=^+ 
(-630 + if 62 + 18 61) 63 + 26 62' + (f 61 - 305) 62 + h^ - 70 61) v^ 
+^h' + (if 62 + 225 + f 61) 63 + if 62^ + (f 61 + iff ) 62 + 25 61 
+1 61')) (cos (v)f + (16 (62 (9 63 + 4 62 + 61)' t;6 _ 14 (9 63 + 4 62 + 61) 
(9 63 + f 62 + 61) v^ + (729 63^ + (-1260 + 162 61 + ifi 62) 63 + 234 62' 
+ (-535 + if 61) 62 - 140 61 + 9 61^) t;2 - 350 - 15 61 - if 62 - 135 63) 
(5 + (963 + 462 + 61) v^) sin {v) + 288t; ((61 + f 63 + fl 62) 
(963 + 462 + h{fv^ - 5 (963 + 462 + 61) (963^ + (61 + 762 - f ) 63 
+1 62' + (i 61 - f) 62 - 2 61) v^ + 225 63' + (25 61 + if^ + 175 62) 63 
+ if 62^ + (f 61 + iff) 62 + 25 61)) cos(t;) + ((61 - 2763) 
(963 + 462 + 6i)%6 4. 63 (f 63 + 61 + ff 62) (963 + 462 + 61) v^ 
+ (-9720 63^ + (4095 - 7560 62 - 1080 61) 63 - 1440 62^ 
+ (-360 61 + 2120) 62 + 555 61) v"^ + 1800 63 + 600 62 + 1325) 
(5 + (963 + 462 + 61) v^) sin iv) + 48 1; ((f 62 + 963 + 61) 
(963 + 462 + 6i)%4 + i (963 + 462 + 61) 

(81 63' + (f 60 + 180 + 27 62 + 18 61) 63-4 62' + (3 61 + 130 + 10 60) 62 
+ (61 + I 60 + 20) 61) t;2 - if 63' + (-if 62 - ^ 60 + 225 - 45 61) 63 
+10 62' + (-25 60 + if - f 61) 62 - I (-10 + I 60 + 61) 61)) 
/(5 + (9 63 + 4 62 + 6i)t;2)' 

After solving the system: 

PL = 0, PL' = 0, PL" = 0, PL'" = 
we get the coefficients: 
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bl = ^ , 02- h - ' 



3D ' 4D ' 2D ' ^ 12D ' (12) 

where D = v^ (cos (f ) + 1) (sin (v)) and 

&o,n«m = 192 (cos(t;))^t;^ — 126 sin(-u)(cos(-u))^t;^ + 99 sin(-u)(cos(-u))^w 

— 126 sin(f)(cos(f ))^f^ — 18 sin(f )f •^ cos(f ) + 630 sin(f )(cos(f))'^w 

— 144 sin(v)(cos(v))^v + 48 sin(w)v^(cos(w))'^ — 288 sin(w)(cos(f ))^f 
-144 (cos(t;))2 - 12 (cos(t;))3 + 336 (cos(t;))^ + 48 sin(w)t;3(cos(w))6 
+30t;2 + 36 cos(t;) + 144 {cos{v)yv^ - 24 (cos(t;))5 + 249 (cos(u))V 
-418 (cos(t;))3t;2 - 662 (cos(t;))S2 + 148 (cos(t;))5t;2 - 99 cos(v) sin(t;)i; 
—9 sin(t>)t> + 66 cos(f)f^ + 96 sin(t>)t>^(cos(f ))^ + 96 sin(t>)(cos(f ))'*t;^ 
-126 sin(t;)(cos(v))S - 18 sin(t;)t;3 - 192 (cos(t;))^ + 176t;2(cos(t;))^ 
-288 sin(-u)(cos(v))^t; 

&i,n«m = 12 + 352 (cos(t;))^t;2 - 36 sin(t;)(cos(t;))V + 168 sin(t;)(cos(w))2i; 
—100 sin(t>)(cos(f))^f^ — 92 sin(t>)f ^ cos(f ) + 96 sin(f )(cos(f))^f 
-672 sin(t;)(cos(t;))S - 384 {cos{v)y + 36 (cos(t;))2 
-48 (cos(t;))3 - 48 (cos(t;))^ + 128 sin(t;)t;3(cos(t;))6 + 33^;^ - 48 cos(t;) 
+448 (cos(t;)) V + 480 (cos(t;))5 - 129 {coslv))^^ - 196 (cos(^;)) V 
—316 (cos(f ))'^f ^ — 464 (cos(f ))^f ^ + 12 cos(f ) sin(f )f — 45 sin(w)t; 
+197 cos(t;)t;^ + 128 sin(t;)t;3(cos(t;))^ - 32 sin(t;)(cos(t;))V 
+504 sin(t>)(cos(f ))S + 4 sin(f )f^ — 288 sin(t>)(cos(f ))^f 

b2,num = 152 (cos(i;))S2 _ 96 (cos(t;))6 + 48 sin(t;)t;3(cos(t;))^ 

-192 sin(w)(cos(v))^t; + 104 (cos(t;))^t;^ - 72 sin(t;)(cos(t;))S 

+48 sin(t;)(cos(t;))V - 244 (cos(t;)) V + 144 (cos(t;))^ 

+216 sin(v)(cos(w))^t; — 44 sin(t;)(cos(t;))^t;^ — 134 (cos(f ))^f ^ 

— 12 (cos(v))^ + 39 sin(f)(cos(w))^f — 44 sin(f)(cos(v))^f^ 

—48 (cos(t;))^ + 79 (cos(f ))^f ^ — 33 cos(f ) sin(f )f — 4 sin(t>)t>^ cos(f ) 

+ 18 cos(f)f^ + 12 cos(f) — 3 sin(f)f — 4 sin(f)f^ + 10 f^ 

b3,num = 208 (cos(t;))^t;2 - 96 (cos(t;))^ - 216 sin(t;)(cos(t;))S 
+120 (cos(t;))S2 + 96 sin(t;)(cos(w))V - 360 (cos(t;))V 
+ 144 (cos(t;))^ - 72 sin(t;)(cos(t;))^t; + 72 sin(w)(cos(t;))V 
+252 sin(t;)(cos(t;))2t; - 157 (cos(t;))2t;2 - 120 sin(^;)(cos(t;))2w3 
— 12 (cos(f ))^ + 149 cos(f )f ^ — 48 cos(f ) + 36 cos(f ) sin(t>)f 
—72 sin(f)f^cos(f) — 45 sin(f)f + 25f^ + 24 sin(f)f^ + 12 
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The Taylor series expansions of the coefficients are given below: 



12629 45767 ^ 27865393 . 557684327 ^ 



3024 9072 11975040 817296480 

235111157089 o 575696865983 ^ 73845973877087 .^ 

v^ + — — v^^ v^^ + ... 



1569209241600 26676557107200 32750603956224000 

20483 45767 . 3549253 , 36881797 . 

bi = v^ H V* v^ 

4032 12096 2280960 99066240 

95714204623 § 138581370311 ^q 106905916402097 ^^ 
^2092278988800^ ~ 35568742809600 ^ ^ 567677135241216000^ 

3937 45767 ^ 3156581 . 21796097 ^ 
^ — ' 'V V H V 



2016 30240 7983360 681080400 
2365857293 § 102137141 ^q 3198002983423 ^^ 



1046139494400 ^ 17784371404800 ^ 283838567620608000 ^ 



17671 45767 ^ 135959 . 14453093 

'J3 



5„ = 11:11^ _ ""'"' ^;2 + ^4 ^6 



12096 181440 47900160 16345929600 

90901339 o 1564247467 ^ 3513993676211 ,0 



896690995200 106706228428800 1703031405723648000 

It is noteworthy that the Taylor series expansions of all four optimized 
methods coincide in the constant term and the coefficient of v'^ and differ on 
the coefficients of f ^ and for higher powers. 

3.4- Error analysis 

We present the principal term of the local truncation error of the five 
methods: 

Classical method: 

PTTF _i^!^„(io)UO 

-r-LJ tjciassical — -^^-„„ V^ ' H 

725760 
Phase ffited method: 

PLTEp,^,,_,,uea = ^^ {/''^ + y^'^co') h'' 

Zero PL and PL' method: 

PLTEu^aer^. = ^^ {v^''^ + UJ^^ + 2 u^'^) h^' 
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Zero PL, PL' and PL" method: 
Zero PL, PL', PL" and PL'" method: 



'; ^ -^ ; 



PLT^3...e„. = ^^ (6i/(^)c.^ + y(^°) + 4yWc.« + 4.;\(«) + c.«y(^)) /. 
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where u is the dominant frequency of the problem. We also present the 
principal term of the local truncation error of the above methods for the case 
of the one-dimensional time-independent Schrodinger equation: 
Classical method: 

PLTEciassicai = -^^ ^'^ \ " 45767^/^5 + 2288351/^4 

7257oU L 
+ ((-2288350 W" - 457670 {W)'^)y - 915340 {W')y')E^ 

+ ((457670 {Wf + 3935962 W^^'^ + 6865050 VT" + 4576700 {W'f)y 

+3661360 {W^^'>)y' + 2746020 {W')y'W)E^ + ((-228835 {WY 

-6865050 IwfW" + (-7871924 W^^'^ - 9153400 {W'f)W - 1327243 W^*^^ 

-9656837 {Wf - 15469246 {W')W^^'>)y - 14645440 {W'){W")y' 

-7322720 W{W^^'^)y' - 2837554 {W^^^)y' - 2746020 {Wf{W')y')E 

+ (45767 {Wf + 2288350 {W)^W" + (3935962 W^W + 4576700 {W'f){Wf 

+ (1327243 W^(6) + 9656837 {Wf + 15469246 {W')W^^^)W 

+2929088 {W')W^^^ + 457671^(8) + 4485166 {W")W^^^ 

+10617944 {W'fW" + 2562952 {W^^^f)y + 915340 {W)%W')y' 

+3661360 {W)\W^^^)y' + (2837554 {W^^^)y' + 14645440 {W'){W")y')W 

+3661360 (W^OV + 366136 (W^(^))i/' + 12814760 (W^")(W^(3))^/ 

+8238060 {W'){W^^'^)y'^ 
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Phase fitted method: 



.10 



-45767W + 45767 W)yE^ 



P LT E Phase-Fitted — _^^_„„ h 

7257oU 

+ ((-1281476 W" + 183068 WW - 183068 {Wf)y - 366136 W'y')E^ 

+ (((4851302 W - 1006874 W)W^" - 274602 {WfW + 274602 {Wf 

+3203690 VrW + 3295224 (W'f)y + ((-549204 W + 1647612 W)W' 

+2562952 W^'^^)y')E'^ + ((-8970332 {W"f + (2013748 WW 

-5858176 (Py)2)W^" + (-7871924 W^+ 1281476 TF)(iy')^ _ 

-14279304 W^W(3) ^ 183068 {WfW - 183068 {Wf + 732272 W^'^W 

-7139652 Vri^(^) - 1281476 Vr(6))i/ - 12448624 VTWV 

+ ((1098408 WW - 2196816 {Wf)W' + 1098408 W^^W 

-6224312 Vri^(3) _ 2562952 W^^^)y')E + ((9656837 Py 

-686505 W) {W"f + (-1006874 (l^ ) W + 4485166 W^ ^ 

+2288350 (HO^ + 10617944 {W'f)W" + (4576700 {Wf_ 

-1281476 VrW)(l^')^ + (29290881^(5) - 1189942^^(3)^7 

+ 15469246 WW^^^)W' + 45767 W^(s) + 45767 {Wf - 45767 {WYW 

+3935962 {WfW^'^^ + (1327243 W^^'^ - 732272 W^^W)W 

+2562952 [W^^^f - Ah767WW^'^^)y + ((-2196816 W 

+14645440 W)W' + 12814760 W^'^^)y'W" + (3661360 {W'f 

+ (8238060 Vr(4) + 915340 {Wf - 549204 {WfW)W' 

+3661360 (W^)W(3) + (2837554 W^(5) - 1098408 W^(3)|F)W^ 

-274602 Wvr(5) + 366136 Vr(^))i/' 



17 



Zero PL and PL' method: 



-'- -10 



PLTEist ,e„. = j:^^;^ h'"" [((-45767 {Wf - 594971 W" - A^mw'' 

+91534 WW)y - 91534 W'y')E^ + ((137301 {Wf - 274602 {WfW 

+ (3157923 W" + 137301 W^)l^ + 2517185 W^^'^ - 1373010 W"W 

+2196816 {W'f)y + (1647612 W^^^ - 549204 W'W + 823806 1^ W')y')E'^ 

+ ((-1235709 W^(6) - 137301 {Wf + 274602 (iy)3W + (-4851302 W" 

-137301 W^)(W^)2 + (-6590448 {W'f + 3386758 W^'W 

-6407380 1^(^))1^ - 320369 ly'W^ - 8283827 (VT")^ + 1373010 WVT^^) 

+2196816 {W'fW - 13089362 W'W^^'^)y + (-2288350 W^^'^ 

-1647612 {WfW + (-51259_04 W^^^) + 1647612 W'W)W - 274602 W^W^ 

-10251808 W'W" + 1830680 WVT^^'))!/')^ + (2929088 W'W^^'^ 

+45767 iy(s) + (1327243 W - 91534 W)W^(6) + 45767 {Wf 

-91534 {WfW + (45767 TF^ + 2288350 W"){Wf + (4576700 {W'f 

-2013748 W"W + 3935962 W^'^'^fWf + (-2562952 (l^')^W^ 

+9656837 {W"f + 320369 VT'W^ - 1464544 WVT^^) 

+15469246 W'W^^^)W + (45767 W^ + 4485166 W")W^^^ 

+10617944 {W'fW" - 1373010 (W^")W+ 183068 {W'fw"^ 

+2562952 {W^'^^f - 2379884 W'W^^W)y + ((-549204 W _ 

+2837554 W)W^^^ + 366136 W^'^^ + 915340 l^'(Vr)3 + (-1098408 W'W 

+3661360 W^(3))(^.)2 ^ (-2196816 WW^(3) ^ 14645440 W^W 

+274602 ly W^)W^ + 8238060 W'W^^^ + (12814760 W" + 183068 W^)W^(3) 
-4393632 W^WW + 3661360 {W'f)y' 



Zero PL, PL' and PL" method: 



1830Q8 yW"E^ + ((45767 (W^)' 



72 X 



j.ndd^tv 725760 
-137301 1^(W^)2 + (17849131^" + 137301 TF^)W^ - 1235709 TFW^" 
-45767 W'^ + 1281476 {W'f + 1876447 W^^)^ ^ 274602 W'y'W 
+ (915340 1^(3) _ 274602 W'W)y')E'^ + ((-91534 {Wf + 274602 {WfW 
+ (-274602 TF^ - 3844428 W^')(W")^ + (-5675108 W^^ + 91534 IF^ 
-5308972 {W'f + 4119030 TFW^")!^ - 11899420 W'W^"^^ - 7597322 {W 
+2746020 {W'fW - 823806 WV" + 1922214 W^^W - 1189942 W^^'^)y 
-1098408 {WfW'y' + (-4027496 W^^^ + 1647612 W'W)y'W 
+ (2196816 W^^W - 8054992 WW" - 549204 W'w'^ - 2013748 W^^'^)y')E 



ll\2 



■2 



+ (45767 {Wf - 137301 (WfW + (137301 W + 2288350 W")(W) 



3 



■3 



+ (-3020622 W W" + 3935962 PyW _ 45757 ^-^ + 4575700 {W'f){W) 



2 



72, 



+ (9656837 (VT")^ + 1327243 VT^^) - 3844428 {W'fW + 961107 VT VT" 
-2196816 W^^W + 15469246 VTW^^'))!^ + 45767 W^^^ - 137301 W^'^W 
+2929088 W^W(5) ^ (4485166 W" + 137301 IF^)W^W + 2562952 (W^(3))2 
-3569826 VTW^^'^W - 2059515 {W'fW + (-45767 TF^ 
+ 10617944 {W'f)W" + 549204 (l^')^W^^)l/ + 915340 {WfW'y' 
+ (-1647612 WW + 3661360 W^'^y)y'{Wf + (14645440 W^W 
-3295224 W^^W + 823806 W^W^ 2837554 W^^^)y'W 
+366136 W^'^'^y' + (-823806 W^^W + 8238060 VTW^^) 

+ (549204 W^ + 12814760 FT") VT^^^ - 91534 PFW^ 
-6590448 WWW + 3661360 {W'f)y' 
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Zero PL, PL', PL" and PL'" method: 

PLI hl?.rd April) = h 



^Srdderiv ^25760 



((1281476 W^^'^ + (-732272 W 

+732272 W)W" +"549204 {W'f)y + 366136 W^%')E^ + ((-1144175 1^(*^) 

+ (2379884 W - 4942836 W)W^^'^ - 10709478 W'W^^'^ - 6910817 {W")^ 

+ (4210564 W^W- 1373010 W^ - 2837554 (1^)2)W^" 

-45767 (W)^ + 183068 {W)W - 274602 {W)W^ 

+ (-4027496 {Wy + 183068 W^)Vr + 2929088 {W')W 

-4:5767W^)y + (-1739146 W^^^ + (2196816 W 

-2929088 W^)W^(3) + io98408 W^WW^' - 5858176 WW" - 549204 PFW^ 

-549204 (iy)2W^')l/')^ + (45767 W^(^) + (1327243 W^ - 183068 W)!^^^) 

+2929088 1^W(5) + (3935962 (Wf + 4485166 VT" + 274602 W'^ 

-2929088 W^TF)l^(^) + 2562952 {W^^^f + (-4759768 WW 

+15469246 W^W^')W"(^) + (-2746020 TF+ 9656837 W^)(W^")2 

+ (2288350 {Wf + 1922214 WW"^ - 4027496 {WfW 

+10617944 (Wy - 183068 W^)W" + 45767 (1^)^ 

-183068 (Vr)W + 274602 {W)W'^ + (-183068 W^^ 

+4576700 (Wy) {Wy + (-5125904 {W'^W + 45767 W^)W^ 

+1098408 (1^')^W^\ + (366136 W^'^^ + (-1098408 W 
+28375541^)1^(5) + 8238060 l^'W^^^) + (3661360 {Wf 
-4393632 W^W+ 1098408 IF^ + 12814760 1^")^"^^^ + (-8787264 W^W 
+ 14645440 W^W^')W"" + 1647612 FTWV' - 2196816 {W^WW 
+3661360 {Wf + 915340 {WfW - 366136 WV')l/' 

The principal terms of the local truncation errors presented above are 
collected in respect to the energy E in descending order. As we can easily 
see, the maximum power of E in the error for each case is: 

• E^ for the classical method 

• E"^ for the phase-fitted method 

• E^ for the zero PL and PL' method 

• E^ for the zero PL, PL' and PL" method and 

• ^2 for the zero PL, PL', PL" and PL'" method. 



/, J. j^ , 
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A low maximum power of E is crucial when integrating the Schrodinger 
equation using a high value of energy. 

3.5. Stability analysis 

The stability analysis of the methods concerns the application of the test 
problem y" = —uoy. 

Here we present the characteristic equations of the five methods: 

C.E.class^cal = 1 + A« + 12k (l^STl s^ " 24192)A7 + ^ (-23622 s^ 

+24192)A6 + ^gg (61449 s^ _ 12096)A5 - ^ s^\* 

+ 12096 (61449 s^ - 12096)A3 + ^ (-23622 s^ + 24192)A2 

+Y2L (17671 s2-24192)A 



C.E.phase-Frtted = "f^ (1 + A^ - 2 A COs(s)) (-|| (A - l)6(cOs(s)) = 



+ 
+ 
+ 
+ 
+ 
+ 



_96_ _|_ _96_ x6 _|_ /'„2 _ 416\ \5 i /'_808 2 , 880\ ^4 , / 1202 2 _ 1120 \ aS 

109 "•" 109 ^ "•" v-^ 109''^ "•" V 327'^ "'"109''^ """ v 327 "^ 109''^ 

808 „2 I 880NX2 I /„2 416n xN/„„„/„N N2 , / 96 96 aG 

~327^ +109''^ +^^ - lOgMJlCOSl^SJj +1-109^109^ 
404 2 _|_ 296n \b 1 /'_480 , 736 2\ \4 , /560 _ 1144 2\ \3 
327 -^ "^109''^ "^ V 109 "^ 327 "^ J^ "T I igg 327 "^ /^ 

480 I 736 „2\ \2 , / 404 2 , 296\ \\ pf.„/'„^ , 32 , 32 ^6 

~T09 + 327^''^ +(>~327^ + T09MJ COS(^SJ + Yog + 109 ^ 

72_ _|_ 137 2\ \5 _|_ (3G_ _ b&_ 2\ \4 , (_30_ 1 302 2\ \3 

109 "•" 327 -^ ;^ "T 1 109 109 ^ /"^ ^ \ 109 """ 327 "^ '"^ 



C.E.,,tDer^. = W (-^ As(A - l)\cos{s)f + ^ (4 A sin(s) 
+s(A2 + 4 A + 1))(A - l)4(cos(s))4 + (-if A (A - If sin(s) 

9 / 48 I 48 \6 192 \ 5 I /'396 , „2n \4 , / 38 „2 504\ \3 
~^ VT25 "•" 125 ^ ^125^ + V125 "•" -^ /^ "'"v 25-^ ~ ^/^ 

+ (fi + s')\' - if A)s)(cos(.))3 + A (-f (A - 1)^ sin(s) + {{s' - i|)A^ 
+ (-§ s^ + 1)A3 + (i .^ - ^)A^ + (-1 s^ + 1)A + .^ 
-i|).)(cos(.))2 + (| A (A - 1)^ sin(s) + if .(§ + | A^ + (-f + .s^)\^ 
(I .^ + f )A^ + (-1 - i s2)A3 + (i .2 + |)A2 + (-f 4 
.SA(A-irsin(s)-l|(l + A^)(| + iA^+(s^-t 



+(| s^ + f )A4 + (-1 - I s')A^ + (| 5' + |)A' + (-f + ^')A)) cos(s) 
-SA(A-irsin(s)-l|(l + A^)(| + iA^+(s^-t)A3 

+ (-tI ^' + W)^' + (^' - W)^)^)(l + A^ - 2 A cos(s))s-i 
(cos(s) + l)"Hcos(s) - 1)"^ 
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C.E.2ndDer^. = -f (-f \H>^ " !)'(«' " 3)(cOs(s))^ 

+f (A^s^ _ I A2 + 3 A sin(s)s + A s^ - 3 A + s^ _ |)A (A - l)2(cos(s))6 
-| (10 s(A2 + 1 + I A)A sin(s) + \^s^ + (4 s^ _ 6)A3 + (18 - 6 s^)\^ 
+(4 s2 - 6)A + s2)(A - l)2(cos(s))5 + (16 (A^ - | A + l)sA (A - 1)^ sin(s) 
+1 AS2 + (f - 16 s2)A5 + (f s2 + f )A^ + (-if -fs' + A s^)\^ 

+ (| s2 ^ ^)^2 + (^ - 16 s2)A + I s2)(cos(s))4 

-4 (-f (A2 + I A + l)sA sin(s) - | AV + (| - f s^)\' + (s^ - ^)\^ 

+(i - f ^^)A - I ^^)(A - l)^(cos(s))3 + i-f (A^ - ^ A + l)sA 

(A - 1)2 sin(s) - f AS2 + (I s2 - f + s4)A5 + (_8 s^ - 16 s^ _ f)\'^ 



-f s'){cos{s)r + 2 (-1 (A^ - I A + l)sA sin(s) - | 



+ (-f s2 + 3 + s^)X^ + (-1 s2 + |)A2 + (-f s2 + 3 + s4)A 
-| s2)(A - 1)2 cos(s) + f s(A2 - i A + 1)A (A - 1)^ sin(s) 
+(1 + A2)(| AS2 + (s^ - I - f s2)A3 + (f + f s2)A2 
+(s4 - I - f s')A + I s'))(l + A2 - 2 A cos(s))s-2 
(cos(s) + l)"2(cos(s) - 1)"^ 

C.E.SrdDer^. = "(1 + A^ - 2 A COs(s))((-f S^ + 32)\%COs{s)y 

-8 (A s(s2 - 6) sin(s) + (4 - f s^)\^ + 3\s^ + A-f .s^)\^ 

{cos{s)f + 12 (s((-5 + s^)\^ + (2 - I s')A - 5 + s^)\ sin(s) 

+(-f s^ + |)A^ + I \h' + ils'- |)A2 + IXs'-fs' + l)X (cos(s))5 

-6 (s((s2 - 3)A4 + (-2 s2 + 4)A3 + (| s' + 2)A2 + (-2 s^ + 4)A 

-3 + s^) sin(s) + I A^s^ + (-8 + f s^)X^ + ("f " ^ ^')^' 

+(-8 + f s2)A + I s2)A (cos(s))4 + (s(AV + (6 - 6 s2)A5 + (66 - 9 s2)A4 

+ (-4 s2 - 3)A3 + (66 - 9 s2)A2 + (6 - 6 s2)A + s^) sin(s) + (30 s^ - Vl)\^ 

+ (_2|5 ^2 _ 4)^4 + (_8 + 1|5 ^2)^3 + (_2|5 ^2 _ 4)^2 ^ (gg ^2 _ 12)A) 

(cos(s))3 + (s(A6s2 + (_21 + 6 s2)A5 + (-9 s^ + f )A4 + (12 s" - 39)A3 
+ (-9 s2 + f )A2 + (-21 + 6 s2)A + s2) sin(s) + (^ s^ + 1)A5 
+(f s2 _ i6)A4 + (i|^ s2 - 2)A3 + (f s=^ - 16)A2 + (^ s^ ^ l)A)(cos(s))2 
+ (-s(AS2 + (3 - 6 s2)A5 + (3 s2 + 9)A4 + (-12 s" + 3)A3 + (3 s^ + 9)A2 
+ (3-6 s2)A + s") sin(s) + (4 - ^ s^^^s + (| 52 ^ 4)^4 + (8 - l|i s2)a3 

^, .,.,.os(s)-(AV-f ^ 

-s2)s(i + A2)sin(s) + (-l- 
+ (-2 - f s2)A3 + 5 A2s2 + (_i _ II s2)A)s-3(cos(sy+ l)-i(sin(s))- 



+ (f S2 + 4)A2 + (4 - ^ S2)^) ^og(g) _ (^4^2 _ ^ ^3 

+(§ + 2 s2)A2 - f A + s2)s(i + A2) sin(s) + (-1 - § s')A5 + 5 A^^ 



12 

From the characteristic equations we evaluate sq and the interval of pe- 
riodicity [0, Sq]. These are given below: 
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• Sq = 0.754 ([0,0.569]) for the classical method 

• So = 0.803 ([0,0.645]) for the phase-fitted method 

• So = 0.874 ([0, 0.763]) for the zero PL and PL' method 

• So = 1.010 ([0, 1.020]) for the zero PL, PL' and PL" method and 

• So = 1.865 ([0,3.478]) for the zero PL, PL', PL" and PL'" method. 

As we can see, by requiring higher derivatives of the phase-lag to be 
vanished, we increase the interval of periodicity, which is a very important 
property. 

4. Numerical results 

4.I. The problems 

The efficiency of the two newly constructed methods will be measured 
through the integration of two real initial value problems with oscillating 
solutions. 

4.I.I. The Schrodinger equation 

The radial Schrodinger equation is given by: 

^"(^) = (^^^ + ^^^) - ^) ^^^) (13) 

where 2 is the centrifugal potential, V{x) is the potential, E is the energy 
and W{x) = 2 +V{x) is the effective potential. It is valid that lim V{x) = 
and therefore lim W{x) = 0. 

We consider E > and divide [0, 00) into subintervals [oj, bi] so that W{x) 
is a constant with value Wi. After this the problem ( IT3l) can be expressed by 
the approximation 

y'i = {W — E) Hi, whose solution is 

y,{x) = A, exp (^^/w-Ex^ + Bi exp (^-^/w-Ex^, (14) 

Ai, Bi G M. 

We will integrate problem (IT3|) with / = at the interval [0, 15] using the 
well known Woods-Saxon potential 
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^r( \ Uo Uiq x-Xo\ 

y{x) = T~^ + 7^7^' g = exp , where (15) 

Uq = —50, a = 0.6, Xo = 7 and Ui = 

a 

and with boundary condition y{0) = 0. 

The potential V{x) decays more quickly than ^^t , so for large x (asymptotic 

region) the Schrodinger equation flT3l) becomes 

y"{x) = (^-^^ - E^ y{x) (16) 

The last equation has two linearly independent solutions kxji{kx) and 
kxni{kx), where ji and rii are the spherical Bessel and Neumann functions. 
When X — s> oo the solution takes the asymptotic form 

y{x) ^ Akx ji{kx) — B kxni{kx) , . 

^ D[sin{kx-TTl/2)+tan{5i) cos{kx - n 1/2)], ^ ■^ 

where 6i is called scattering phase shift and it is given by the following ex- 
pression: 

tan (5 ) - ^^^'^ '^(^^+i) ~ y(^i+^) SJ^i) ^gx 

y{xi+i)C{x.i)-y{x./)C{xi+iy 

where S{x) = kxji{kx), C{x) = kxni{kx) and Xj < Xj+i and both belong 
to the asymptotic region. Given the energy we approximate the phase shift, 
the accurate value of which is 7r/2 for the above problem. 
We will use three different values for the energy: 

• Ei = 989.701916 

• ^2 = 341.495874 

• ^3 = 163.215341 



As for the frequency w we will use the suggestion of Ixaru and Rizea [11 



^_ V£^, .e 10, 6.5] 

^ VE, X G [6.5, 15] 
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4.1.2. The N- Body Problem 

The N-body problem is the problem that concerns the movement of N 
bodies under Newton's law of gravity. It is expressed by a system of vector 
differential equations 



f =G f 7i^I,f> . , = l,2,..,iV (20) 

where G is the gravitational constant, rrij is the mass of body j and yi is the 
vector of the position of body i. 

It is easy to see that each vector differential equation of (|20|) can be 
analyzed into three simplified differential equations, that express the three 
directions x, y,z. So y] —yi expresses the difference between the coordinates 
of bodies j and i for the corresponding direction, while | yj — yi \ represents 
the distance between bodies i and j. 

The above system of ODEs cannot be solved analytically. Instead we 
produce a highly accurate numerical solution by using a 10-stage implicit 
Runge-Kutta method of Gauss with 20th algebraic order, that is also sym- 
plectic and A-stable. The method can be easily reproduced using simplifying 
assumptions for the order conditions (see |7|]). 

The reference solution is obtained by using the previous method to in- 
tegrate the N-body problem for a specific time-span and for different step- 
lengths. 

In order to find the step-length hopt that gives the best approximation, 
we have to keep in mind that the total error of a numerical method that 
integrates a system of ODEs consists of the error due to the truncation error 
of the method and the roundoff error of all computations. While the global 
truncation error of the method tends to zero, while h decreases, the opposite 
happens to the roundoff, which tends to infinity. 

If Vacc is the analytical solution for a specific time-span of the problem, 
then let e„ = \\y,^^ - yacc\\ and £„ = WVh^^, -Vhjl, where y,^^ is the approx- 
imate solution of y using a step- length /i„. e„ represents the actual error 
of the approximation and e„ is the best known approximation to the actual 
error, being the difference of two approximations with different step-lengths. 
We see that, when /i„ — > hopt =^ ^n -^ ^min and En — > £min- The minimum 
values of the errors e,„m and Emin are positive numbers and depend on the 
software that is used for the integration and the computer system that it 
runs on. We can also see that e„ and e„ have similar behavior around riopt, 
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meaning that they increase and decrease simultaneously. According to these 
we find the step-length h^pt that minimizes En, which is easily calculated for 
every hn- 

In [8| the data for the five outer planet problem is given. This system 
consists of the sun and the five most distant planets of the solar system. 
In Table ^T\ we can see the masses, the initial position components and the 
initial velocity components of the six bodies. Masses are relative to the sun, 
so that the sun has mass 1. In the computations the sun with the four inner 
planets are considered one body, so the mass is larger than one. Distances 
are in astronomical units, time is in earth days and the gravitational constant 
is G = 2.95912208286 • 10"^ 



Planet 


Mass 


Initial Position 


Initial Velocity 


Sun 


1.00000597682 












Jupiter 


0.000954786104043 


-3.5023653 


0.00565429 






-3.8169847 


-0.00412490 






-1.5507963 


-0.00190589 


Saturn 


0.000285583733151 


9.0755314 


0.00168318 






-3.0458353 


0.00483525 






-1.6483708 


0.00192462 


Uranus 


0.0000437273164546 


8.3101420 


0.00354178 






-16.2901086 


0.00137102 






-7.2521278 


0.00055029 


Neptune 


0.0000517759138449 


11.4707666 


0.00288930 






-25.7294829 


0.00114527 






-10.8169456 


0.00039677 


Pluto 


1/(1.3-108) 


-15.5387357 


0.00276725 






-25.2225594 


-0.00170702 






-3.1902382 


-0.00136504 



(21) 



The system of equations (l20l) has been solved for t G [0, 10^], for which 
time-span, the previously mentioned method of Gauss produces a 10.5 deci- 
mal digits solution. 

We have used uj = 0.00145044732989, which is the dominant frequency of 
the problem, as evaluated by the square root of the spectral radius of matrix 
A, if the problem is expressed in the form y" = Ay + B. 
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4-2. The methods 



The classical method developed by Quinlan and Tremaine [18 



The phase-fitted method developed by Panopoulos, Anastassi and 



Simos 16 



• The zero PL and PL' method developed here 

• The zero PL, PL' and PL" developed here 

• The zero PL, PL', PL" and PL'" developed here 

4-3. Comparison 

We are presenting the accuracy of the methods expressed by — log^Q (error 
at the end point) versus the log^^Q (total steps). In Figures [H [2] and [3] we 
are presenting the efficiency of the methods for the Schrodinger equation 
using a value for the energy equal to i) 989.701916, ii) 341.495874 and iii) 
163.215341. Also in FigurelHwe present the efficiency for the N-body problem 
and particularly the five outer planet problem. 

We see that for each successive derivative of the phase-lag nullified, we 
gain in efficiency for both IVPs tested here. 

5. Conclusions 

We have developed three new optimized eight-step symmetric methods 
with zero phase-lag and derivatives. We showed that the more derivatives 
of the phase-lag are vanished, the bigger the interval of periodicity and the 
higher the efficiency of the method. This is the case for both problems tested 
here. Also the local error truncation analysis shows the relation of the error to 
the energy, revealing the importance of nullified phase-lag derivatives when 
integrating the Schrodinger equation, especially when using high value of 
energy. 
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